/*
 * GaussianNoise.cpp
 *
 *  Created on: Jun 10, 2011
 *      Author: nzhao
 */

#include "share/Usrdef.h"
#include "share/parameters.h"
#include "special/spectrum.h"
#include "core/Noise.h"
#include "share/MatAlgebra.h"
#include "share/gaussian_integration.h"

//////////////////////////////////////////////////////////////////////////
//Project Path:: toggle between Local and Server machines

char	CODE_LOCAL[] = "/home/nzhao/research/codes/gaussian-noise/";
char	CODE_SERVER[] = "/home/user/nzhao/Projects/GaussianNoise/";
char	CODE_PATH[200];
char	WORKSHOP_LOCAL[] = "/home/nzhao/research/projects/AntiCrossing/";
char	WORKSHOP_SERVER[] = "/home/user/nzhao/Projects/AntiCrossing/";
char	WORKSHOP[200];
char 	DEBUG_PATH[] = "data/output/dbg/";

//////////////////////////////////////////////////////////////////////////
//Global Variables
clock_t     START, FINISH, ABS_TIME_START, ABS_TIME_END;
double      TIME_USED, ABS_TIME_USED;

noise_ensemble b_t;


//////////////////////////////////////////////////////////////////////////
//Global function declarations
void Run_Time(const char * flag);
void path_samples();
void check_correlation_fucntion();
void test_phase_calculations(int path_idx);
void coherence(double * coh_re, double * coh_im);
void export_results(double * corr_func, double * coh_re, double * coh_im);

int main(int argc, char **argv)
{
	Run_Time("start");
	int i, j;

	if (strcmp(WorkshopPath, "LOCAL") == 0)
	{
		strcpy(WORKSHOP, WORKSHOP_LOCAL);
		strcpy(CODE_PATH, CODE_LOCAL);
	}
	else
	{
		strcpy(WORKSHOP, WORKSHOP_SERVER);
		strcpy(CODE_PATH, CODE_SERVER);
	}
	cout << "Workshop path = " << WORKSHOP << endl;
	cout << "CODE path = " << CODE_PATH << endl;

	ParameterResolve(argc, argv);
	Generate_Power_List();

	b_t.initialize(MIN_TIME, MAX_TIME, NTIME, MIN_FREQ, MAX_FREQ, NFREQ, NPATH, N_INTERPOLATION);
	b_t.generate_bath(SIM_SEED);

    double * corr_func = new double [NTIME];
    b_t.calc_correlation_func(corr_func);

//	test_phase_calculations(66);
//	path_samples();

	double * coh_re, * coh_im;
	coh_re = new double [NTIME];
	coh_im = new double [NTIME];

	coherence(coh_re, coh_im);

	Run_Time("end");

	export_results(corr_func, coh_re, coh_im);
}

void Run_Time(const char * flag)
{
    //if(k==1)    //start
    if(strcmp(flag, "start") == 0 )
    {
        START = clock();
        time(&ABS_TIME_START);
    }
    else if(strcmp(flag, "end") == 0 )
    {
       	FINISH = clock();
       	time(&ABS_TIME_END);
        TIME_USED = (double)(FINISH - START) / CLOCKS_PER_SEC;
        ABS_TIME_USED = difftime(ABS_TIME_END, ABS_TIME_START);
        cout << "#######################################" << endl;
        cout << "calculation time:" << endl;
        cout << "CPU time = " <<  TIME_USED << " second;" << endl;
        cout << "Absolute Time = " << ABS_TIME_USED << " second." <<   endl;
    }
    else
    {
        cout << "do nothing!" << endl;
    }
}

void coherence(double * coh_re, double * coh_im)
{
	int i, n;
	double c0, c1, c2, dz, zval;

	double * phase_mat;
	phase_mat = new double [NTIME * NPATH];

	for(i = 0; i < NPATH; i ++)
	{
		if(i % 100 == 0)
			cout << "calculate phase of " << i << " th path." << endl;
		b_t.calc_phaseFull_gsl_int(phase_mat + i * NTIME, i);
	}

	b_t.set_phase_mat(phase_mat);
	b_t.calc_dephasing(coh_re, coh_im );

	cout << "#######################################" << endl;
	for(i = 0; i < NTIME; i++)
		cout << i*TIME_STEP <<"\t" << coh_re[i] << "\t" << coh_im[i] << "\t" << sqrt(coh_re[i]*coh_re[i] + coh_im[i]*coh_im[i])<< endl;

	delete [] phase_mat;
}

void export_results(double * corr_func, double * coh_re, double * coh_im)
{
	int i;
    char coh_file[200];
    char * war;

    strcpy(coh_file, WORKSHOP);
    strcat(coh_file, "data/output/");
    if(IS_WARNING > 0)
    {
    	war = itoa(IS_WARNING, 10);
    	strcat(coh_file, "WARNING_");
    	strcat(coh_file, war);
    	strcat(coh_file, "_");
    }

    strcat(coh_file, OUTPUT_FILE);

    ofstream coh(coh_file);
    if(!coh) assert(0);


    coh << "#info\t" << "Start Time: " << ctime(&ABS_TIME_START) ;
    coh << "#info\t" << "Finish Time: " << ctime(&ABS_TIME_END) ;
    coh << "#info\t" << "Absolute Time: " << ABS_TIME_USED << " seconds." << endl;
    coh << "#info\t" << "CPU Time: " << TIME_USED << " seconds." << endl;
    coh << "#info\t" << "VERSION: " << VERSION_INFO << endl ;
    if(IS_WARNING > 0)
    	coh << "#info\t" << "MAX_ERROR: " << setiosflags(ios::scientific) << setprecision(10)<< MAX_ERROR << endl;

    coh << "#para\t" << "Start time = " << MIN_TIME << endl ;
    coh << "#para\t" << "End time = " << MAX_TIME << endl ;
    coh << "#para\t" << "Ntime = " << NTIME << endl ;
    coh << "#para\t" << "MIN Freq = " << MIN_FREQ << endl ;
    coh << "#para\t" << "Max. Freq = " << MAX_FREQ << endl ;
    coh << "#para\t" << "NFreq = " << NFREQ<< endl ;
    coh << "#para\t" << "N_Interpolation = " << N_INTERPOLATION<< endl ;
    coh << "#para\t" << "Noise amplitude = " << NOISE_AMPLITUDE << endl ;
    coh << "#para\t" << "Correlation time = " << CORRELATION_TIME << endl ;
    coh << "#para\t" << "Npath = " << NPATH << endl ;
    coh << "#para\t" << "seed = " << SIM_SEED<< endl ;
    coh << "#para\t" << "Omega z = " << OMEGA_Z<< endl ;
    coh << "#para\t" << "Omega v = " << OMEGA_V << endl ;
    coh << "#para\t" << "Abs eps = " << ABS_EPSLON << endl ;
    coh << "#para\t" << "Gauss Interval = " << GAUSS_INTERVAL << endl ;
    coh << "#para\t" << "Pulse Seq = " << STR_SEQ_NAME << endl ;
    coh << "#para\t" << "Pulse Num.= " << PULSE_NUM << endl ;

    coh << "time" << "\t" << "correlation" << "\t" << "re" << "\t" << "im" << "\t" << "abs" << endl;
    coh << setw(12) << setiosflags(ios::scientific) << setprecision(10);
    for(i = 0; i < NTIME; i++)
    {
    	coh << i*TIME_STEP << "\t" << corr_func[i] << "\t" << coh_re[i] << "\t" << coh_im[i] << "\t" << sqrt(coh_re[i]*coh_re[i] + coh_im[i]*coh_im[i])<< endl;
    }

	coh.close();
}

double b_t_val(double tval, void * para)
{
	double res;
	int n = * (int *) para;

	res = b_t.eval(n, tval);
	return res;
}

void path_samples()
{
    char b_file[200];

    strcpy(b_file, WORKSHOP);
    strcat(b_file, DEBUG_PATH);
    strcat(b_file, "Bt_");
    strcat(b_file, PARA_CHAR);
    strcat(b_file, ".dat");

    ofstream btSamples(b_file);
    if(!btSamples) assert(0);

    int n;
    double t, tmin = 0.0, tmax =100.0, dtval = 0.1;
	for(t = tmin; t < tmax; t+= dtval)
	{
		btSamples << t << "\t";
		for(n = 0; n < 10; n ++)
			btSamples << b_t.eval(n, t) << "\t";
		btSamples << endl;
	}
	btSamples.close();
}

void test_phase_calculations(int path_idx)
{
	int i;
	double * phase1_exact, * phase_gsl, * phase_gsl_interpolation;
	double * phase2_exact, * phase_gsl_interpolation2;
	double * phaseFull_gsl_interpolation;

	phase1_exact = new double [NPATH * NTIME];
	phase_gsl= new double [NTIME];
	phase_gsl_interpolation= new double [NTIME];
	phase_gsl_interpolation2 = new double [NTIME];

	phase2_exact = new double [NTIME];
	phaseFull_gsl_interpolation = new double [NTIME];

	b_t.calc_phase1(phase1_exact);
//	b_t.calc_phase1_gsl_int(phase_gsl, path_idx);
	b_t.calc_phase1_gsl_int1(phase_gsl_interpolation, path_idx);

	cout << setw(12) << setiosflags(ios::scientific) << setprecision(10);
	for(i = 0; i < NTIME; i ++)
	{
		cout << i*TIME_STEP << "\t"
				<< phase1_exact[i * NPATH + path_idx] << "\t"
//				<< phase_gsl[i] << "\t"
				<< phase_gsl_interpolation[i] << "\t dif1 = "
//				<< phase_gsl[i] - phase_exact[i * NPATH + path_idx] << "\t diff2 = "
				<< phase_gsl_interpolation[i] - phase1_exact[i * NPATH + path_idx] << endl;
	}

	b_t.calc_phase2(phase2_exact, path_idx);
	b_t.calc_phase2_gsl_int(phase_gsl_interpolation2, path_idx);

	cout << setw(12) << setiosflags(ios::scientific) << setprecision(10);
	for(i = 0; i < NTIME; i ++)
	{
		cout << i*TIME_STEP << "\t"
				<< phase2_exact[i ] << "\t"
				<< phase_gsl_interpolation2[i] << "\t diff = "
				<< phase_gsl_interpolation2[i] - phase2_exact[i] << endl;
	}

	b_t.calc_phaseFull_gsl_int(phaseFull_gsl_interpolation, path_idx);
	cout << setw(12) << setiosflags(ios::scientific) << setprecision(10);
	for(i = 0; i < NTIME; i ++)
	{
		cout << i*TIME_STEP << "\t" << phaseFull_gsl_interpolation[i] << endl;
	}
}
